library(tidyverse)
library(hablar)
library(patchwork)
library(plotly)
library(ggparallel)
library(ggforce)
library(lubridate)
library(fmsb)
library(rstatix)
#functions
sd.prop <- function(x, n){
sqrt((x*(1-x))/n)
}
se <- function(x){
sd(x)/sqrt(14*729) #n=No of years*number of parameter combinations for 2004-2018
}
shap_test <- function(df){
shapiro.test(df$med_DBH)
}
shap_test2 <- function(df){
shapiro.test(df$aver)
}
wil_test <- function(df){
wilcox.test(med_DBH~PFT, df, alternative='less')
}
wil_testG <- function(df){
wilcox.test(aver~PFT, df, alternative='greater')
}
#Not significantly different
wil_testNS <- function(df){
wilcox.test(aver~PFT, df)
}
#dropping year function
drop_y <- function(df){
select(df,-Year)
}
#Kappa function
fun_Kap2 <- function(x){
Kappa.test(x)
}
#getting Kappa values
est_Kap <- function(x){
x$Result$estimate
}
#getting p values
p_Kap <- function(x){
x$Result$p.value
}
#effect size
wil_test_eff <- function(df){
df %>% pivot_longer(cols=c(freq,freqW), names_to='type', values_to='No') %>%
wilcox_effsize (No~type)
}
Analysis and statistical code is shown for RB+Pollen (named as Polim_)
#Field data/Obs Monitoring1000 (70m*150m)
uryu<- read.csv('uryu.csv',row.names = 1,
colClasses = c(Family='factor',
Latin_name='factor',
PFT='factor'))
#Monitoring1000 for 1 hectare, biomass in kg
uryuYear <- read.csv('uryuYear.csv',row.names = 1)
#Seed
#1983-
Yacorns <- read.csv('Yacorns.csv', row.names = 1)
#2004
Yseed <- read.csv('Yseed.csv', row.names = 1,
colClasses = c(Latin_name='factor',
mastingW='factor',
mastingN='factor'))
#Simulation data
#QuercusNo and FloweringNo at the time of flowering, Masting and Masting_n at end of year
#MassRatio and FloweringRatio at time of flowering
Polim <- read.csv('Polim.csv', row.names = 1)
#Field data stats
uryuYear %>%
summarise(meanNo=mean(TreeNo), sdNo=sd(TreeNo),
meanNoProp=mean(TreeProp), seNoProp=mean(sdNoProp),
meanBiom=mean(Biomass), sdBiom=sd(Biomass),
meanBiomProp=mean(BiomassProp), seBiomProp=mean(sdBiomProp))
#Simulation between 2005-2018
#RB+Polim
#tree number
Polim %>% filter (Year>2004) %>%
group_by(Rc, Li, Bu, Year) %>% #will calculate mean from 5 runs for every year
summarise(Tree=mean(TreeNo), Percent=mean(NoRatio)) %>% ungroup() %>%
mutate(sePerc=sd.prop(x=Percent, n=Tree)) %>%
summarise(mean=mean(Tree), se=se(Tree), meanPerc=mean(Percent), sePercM=mean(sePerc))
#Biomass
Polim %>% filter (Year>2004) %>% group_by(Rc, Li, Bu, Year) %>%
summarise(biomass=mean(WoodyPftMass), Quercus=mean(Masting), Percent=Quercus/biomass,
Tree=mean(TreeNo), sePerc=sd.prop(x=Percent, n=Tree)) %>%
ungroup() %>%
summarise(mean=2*mean(biomass), se=2*se(biomass), meanQ=2*mean(Quercus), seQ=2*se(Quercus),
PercentM=mean(Percent), PercentSE=mean(sePerc))
#Obs
uryu %>% group_by(PFT) %>%
summarise(biomassPft=mean(biomass),
DBHmin=range(dbh)[1],
DBHmax=range(dbh)[2],
DBHmed=median(dbh),
DBHav=mean(dbh)) #Masting is 642 kg
#normality of biomass and dbh
uryu2 <- uryu %>%
filter(PFT %in% c('Masting', 'TeBS')) %>%
group_by(Year,PFT) %>% summarise(MeanB=mean(biomass), MedBH=median(dbh))
tapply(uryu2$MeanB,uryu2$PFT, shapiro.test)#both normal
## $BoBS
## NULL
##
## $BoNE
## NULL
##
## $Masting
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.92512, p-value = 0.2604
##
##
## $TeBS
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.98663, p-value = 0.997
tapply(uryu2$MedBH,uryu2$PFT, shapiro.test)#both normal
## $BoBS
## NULL
##
## $BoNE
## NULL
##
## $Masting
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.93776, p-value = 0.3904
##
##
## $TeBS
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.96244, p-value = 0.763
#t-test
t.test(MeanB~PFT, data=uryu2) #mean biomass
##
## Welch Two Sample t-test
##
## data: MeanB by PFT
## t = 40.96, df = 13.091, p-value = 3.297e-15
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 508.4844 565.0672
## sample estimates:
## mean in group Masting mean in group TeBS
## 644.1987 107.4229
t.test(MedBH~PFT, data=uryu2) #median DBH
##
## Welch Two Sample t-test
##
## data: MedBH by PFT
## t = 5.7635, df = 21.297, p-value = 9.6e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.674635 1.435305
## sample estimates:
## mean in group Masting mean in group TeBS
## 11.11015 10.05518
#Sim
#RB+Pollen-median DBH
trial2 <- Polim %>% select(Rc:run, TeNE_medDBH:BoBS_medDBH) %>%
pivot_longer(cols=TeNE_medDBH:BoBS_medDBH, names_to='PFT', values_to='median') %>%
group_by(Year, PFT, Rc, Li, Bu) %>%
summarise(med_DBH=mean(median)) #eliminating runs
#normality
#nesting
trial2 %>% mutate(med_DBH=med_DBH+rnorm(n(),mean = 0, sd=0.0001))%>% #correcting for all equal problem
group_by(PFT, Rc, Li, Bu) %>%
nest() %>%
mutate(shapTest=map(data, shap_test), glance=map(shapTest, broom::glance)) %>%
unnest(glance) %>% filter(p.value<0.05) %>%
group_by(PFT) %>% count() #Mostly non-normal
#Wilcox test
trial_Wil <- Polim %>% select(Rc:run, TeNE_medDBH:BoBS_medDBH) %>%
pivot_longer(cols=TeNE_medDBH:BoBS_medDBH, names_to='PFT', values_to='median') %>%
filter(PFT %in% c('TeBS_medDBH','Masting_medDBH')) %>%
mutate(PFT=fct_relevel(PFT, 'TeBS_medDBH')) %>%
group_by(Year, PFT, Rc, Li, Bu) %>% summarise(med_DBH=mean(median))
trial_Wil %>% group_by(Rc, Li, Bu) %>%
nest() %>%
mutate(WilTest=map(data, wil_test), glance=map(WilTest, broom::glance)) %>%
unnest(glance) %>% filter(p.value < 0.05) #Masting PFT never higher
#RB+Pollen Mean Biomass
Polim_stat <- Polim %>% mutate(TeNE=TeNE/TeNE_n,
TeBS=TeBS/TeBS_n,
BoNE=BoNE/BoNE_n,
BoNS=BoNS/BoNS_n,
BoBS=BoBS/BoBS_n,
Masting=Masting/Masting_n) %>%
select(Rc:run, TeNE:Masting) %>%
pivot_longer(cols=TeNE:Masting, names_to='PFT', values_to='aver') %>%
filter(PFT %in% c('TeBS','Masting'))
#nesting
Polim_stat %>% group_by(PFT, Rc, Li, Bu) %>% #taking together 5 runs Masting or TeBS
nest() %>%
mutate(shapTest=map(data, shap_test2), glance=map(shapTest, broom::glance)) %>%
unnest(glance) %>% filter(p.value<0.05) %>%
group_by(PFT) %>% count() #Mostly non-normal
#Getting parameter combinations, Masting S >
Polim_par_g <- Polim_stat %>% group_by(Rc, Li, Bu, Year, PFT) %>% summarise(aver=mean(aver)) %>%
group_by(Rc, Li, Bu) %>%
nest() %>%
mutate(WilTest=map(data, wil_testG), glance=map(WilTest, broom::glance)) %>%
unnest(glance) %>% filter(p.value<0.05) %>%
unite('parameter', c('Rc', 'Li','Bu'), sep='_', remove = F) %>% .$parameter
#NS
Polim_par_e <- Polim_stat %>% group_by(Rc, Li, Bu, Year, PFT) %>% summarise(aver=mean(aver)) %>%
group_by(Rc, Li, Bu) %>%
nest() %>%
mutate(WilTest=map(data, wil_testNS), glance=map(WilTest, broom::glance)) %>%
unnest(glance) %>% filter(p.value>=0.05) %>%
unite('parameter', c('Rc', 'Li','Bu'), sep='_', remove = F) %>% .$parameter
#parameters unite
Polim_par <- union(Polim_par_e, Polim_par_g)
#filtered dataset
Polim_f <- Polim %>%
unite('parameter', c('Rc', 'Li','Bu'), sep='_', remove = F) %>%
filter(parameter %in% Polim_par)
Polim_f %>% filter(Year>2004) %>%
group_by(Rc, Li, Bu, run) %>% summarise(FlowerCV=sd(FloweringRatio)/mean(FloweringRatio)) %>%
group_by(Rc, Li, Bu) %>% summarise(meanFlowerCV=mean(FlowerCV), seFlowerCV=sd(FlowerCV)/sqrt(5)) %>%
arrange(desc(meanFlowerCV))
#Uryu
e <- uryu %>% group_by(Year, PFT) %>%
summarise(biomass=sum(biomass), count=n()) %>% #for every Year, every family average
mutate(aver=biomass/count) %>%
ggplot(aes(x=PFT, y=aver))+
geom_boxplot()+
labs(title='(e)', y='Mean biomass/tree (kg dry mass)')+
theme_bw()
#RB+Pollen
c <- Polim %>% mutate(TeNE=TeNE/TeNE_n,
TeBS=TeBS/TeBS_n,
BoNE=BoNE/BoNE_n,
BoNS=BoNS/BoNS_n,
BoBS=BoBS/BoBS_n,
Masting=Masting/Masting_n) %>%
select(Rc:run, TeNE:Masting) %>%
pivot_longer(cols=TeNE:Masting, names_to='PFT', values_to='aver') %>%
filter(Year>2004) %>% group_by(Year, PFT, Rc, Li, Bu) %>%
summarise(aver=mean(aver)) %>%
ggplot(aes(x=PFT, y=aver*1000*2))+
geom_boxplot()+
labs(title='(c)', y='Mean biomass/tree (kg dry mass)')+
theme_bw()
c+e
#RB+Pollen
Polim_mast_meanBiom <- Polim %>% mutate(TeNE=TeNE/TeNE_n,
TeBS=TeBS/TeBS_n,
BoNE=BoNE/BoNE_n,
BoNS=BoNS/BoNS_n,
BoBS=BoBS/BoBS_n,
Masting=Masting/Masting_n) %>%
select(Rc:run, TeNE:Masting) %>%
pivot_longer(cols=TeNE:Masting, names_to='PFT', values_to='aver') %>%
filter(Year>2004) %>% group_by(Year, PFT, Rc, Li, Bu) %>%
summarise(aver=mean(aver)) %>% filter(PFT=='Masting')
trial <- Polim_mast_meanBiom %>% filter(Rc==5, Year==2010) %>% ungroup() %>%
select(Bu, Li, aver) %>% mutate(aver=aver*1000*2) %>%
convert(num(Bu, Li)) %>% mutate(Bu=Bu*100, Li=Li*100)
trial2 <- xtabs(aver~Li+Bu, data=trial) #9*9
trial3 <- as.numeric(trial2)
trial4 <- matrix(trial3, nrow = 9, ncol=9)
plot_ly(x=as.numeric(colnames(trial2)),
y=as.numeric(rownames(trial2)),
z=trial4, type = 'surface', colorscale = list(c(0, 1), c('white', "blue")),
cmin=0, cmax=220,
contours=list(
z=list(show = TRUE, start = 0.0, end = 220, size = 6, color = 'black')
)) %>%
layout(title=list(text=paste0('(c)'), y=0.97, x=0.2),font=list(size=14),
scene=list(aspectratio=list(x=1,y=1,z=1),
xaxis=list(title='<i>B<sub>u</sub></i> (%)'),
yaxis=list(title='<i>L<sub>T</sub></i> (%)'),
zaxis=list(title='Mean biomass (kg)',range=c(0,220)),
camera = list(eye = list(x = 1.3, y = 1.3, z = 1.3),
center = list(x=0, y=0, z=-0.2)))) %>% hide_colorbar()
#other filtered simulations
Basic_f <- read.csv('Basic_f.csv',row.names = 1)
BasicW_f <- read.csv('BasicW_f.csv',row.names = 1)
PolimW_f <- read.csv('PolimW_f.csv',row.names = 1)
#parameters filtered for each masting function
Basic_f <- Basic_f %>% mutate(Type='RB')
BasicW_f <- BasicW_f %>% mutate(Type='RB-W')
Polim_f <- Polim_f %>% mutate(Type='RBP')
PolimW_f <- PolimW_f %>% mutate(Type='RBP-W')
#Combining
combined <- bind_rows(Basic_f, BasicW_f, Polim_f, PolimW_f)
combined <- combined %>% select(Rc:Bu, Type) %>% convert(fct(Type)) %>%
mutate(Type=fct_relevel(Type, 'RB', 'RBP', 'RB-W'))
combined2 <- combined %>% convert(num(Rc:Bu)) %>%
mutate(Li=Li*100, Bu=Bu*100) %>%
group_by(Rc, Li, Bu, Type) %>% count()
combined3 <- gather_set_data(combined2,1:3)
combined3 %>% mutate(Type=fct_recode(Type, '(a)'='RB',
'(b)'='RB-W',
'(c)'='RBP',
'(d)'='RBP-W')) %>%
ggplot(aes(x, id = id, split = y, value = n)) +
geom_parallel_sets( alpha = 0.3, axis.width = 0.1) +
geom_parallel_sets_axes(axis.width = 0.3, fill='white', color='grey') +
geom_parallel_sets_labels(colour = 'black' , angle = 0)+
facet_wrap(~Type)+
theme_minimal()+
theme(axis.text.y=element_blank(),
axis.ticks = element_blank(),
axis.text.x = element_text(size=13, face='bold'),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
strip.text.x = element_text(size=12, hjust = 0.1))+
scale_x_discrete (name='Parameters',
labels = c(expression(paste(italic('B'[u]),' (%)')),
expression(paste(italic('L'[T]),' (%)')),
expression(paste(italic('R'[c])))))
# RB+Pollen
trialPolBu3 <- Polim_f %>% filter(Rc==3, Bu==0.03, Year>1983) %>%
group_by(Year, Li) %>% summarise(FloweringRatio=mean(FloweringRatio)) %>%
ungroup() %>%
select(Year, Li, FloweringRatio) %>%
convert(num(Li)) %>% mutate(Li=Li*100)
trialPolBu4 <- Polim_f %>% filter(Rc==3, Bu==0.04, Year>1983) %>%
group_by(Year, Li) %>% summarise(FloweringRatio=mean(FloweringRatio)) %>%
ungroup() %>%
select(Year, Li, FloweringRatio) %>%
convert(num(Li)) %>% mutate(Li=Li*100)
plot_ly(data=trialPolBu3[which(trialPolBu3$Year>2004),], x = ~Li, y = ~Year, z = ~FloweringRatio, split = ~Li, type = 'scatter3d',
mode = 'lines',
line = list(color='black',width = 2)) %>% #works really well, different color for lines
add_trace(data = trialPolBu4[which(trialPolBu4$Year>2004),],
x = ~Li, y = ~Year, z = ~FloweringRatio, split = ~Li,
line = list(color='blue',width = 2)) %>%
layout(font=list(size=15), title=list(text='(c)',y=0.9, x=0.2),
scene=list(aspectratio=list(x=1.2,y=1,z=1),
xaxis=list(title='<i>L<sub>T</sub></i> (%)', autotick=F,
dtick=2.5, range=c(2, 23)),
yaxis=list(title='Year'),
zaxis=list(title='Flowering synchrony', autotick=F,
dtick=0.1, range=c(0.01,0.58)),
camera = list(eye = list(x = 1.3, y = 1.3, z = 1.3),
center = list(x=0, y=0, z=-0.25)))) %>% hide_legend()
Polim_f %>% filter(Year>1983) %>%
group_by(Rc, Li, Bu, run) %>% summarise(CV=sd(AcornW)/mean(AcornW)) %>%
group_by(Rc, Li, Bu) %>% summarise(meanCV=mean(CV), seCV=sd(CV)/sqrt(5)) %>%
arrange(desc(meanCV)) %>% head(5) %>% left_join(Polim_f) %>%
group_by(Year, parameter, meanCV) %>%
summarise(Acorn=mean(AcornW), SD_Acorn=sd(AcornW)) %>%
convert(fct(parameter)) %>%
mutate(parameter=fct_reorder(parameter, desc(meanCV))) %>%
filter(Year>1983) %>%
group_by(parameter) %>%
summarize(ratio=SD_Acorn/Acorn, Acorn=(Acorn-mean(Acorn))/sd(Acorn), #ratio to scale later SD
SD_Acorn=(SD_Acorn-mean(SD_Acorn))/sd(SD_Acorn), #standardized to SD
SD_Acorn=SD_Acorn*ratio, Year=Year) %>%
ggplot(aes(x=Year))+
geom_ribbon(aes(fill=parameter,
ymin=Acorn-SD_Acorn,
ymax=Acorn+SD_Acorn), alpha=0.1)+
geom_line(aes(y=Acorn, color=parameter))+
labs(title='(c)', y='Acorn volume (SD)')+
theme_bw()+
theme(legend.title = element_blank(),
legend.position = c(0.5, 0.88),
legend.background=element_rect(fill='transparent', size=NULL),
legend.text.align = 0,
axis.title.y = element_blank())+
scale_fill_discrete(labels=c(
'3_0.175_0.04'=
expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=17.5%,',italic('B'[u]),'=4% CV=0.65')),
'3_0.05_0.04'=
expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=5%,',italic('B'[u]),'=4% CV=0.64')),
'3_0.025_0.04'=
expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=2.5%,',italic('B'[u]),'=4% CV=0.63')),
'3_0.2_0.04'=
expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=20%,',italic('B'[u]),'=4% CV=0.62')),
'3_0.15_0.04'=
expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=15%,',italic('B'[u]),'=4% CV=0.61'))))+
scale_color_discrete(labels=c(
'3_0.175_0.04'=
expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=17.5%,',italic('B'[u]),'=4% CV=0.65')),
'3_0.05_0.04'=
expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=5%,',italic('B'[u]),'=4% CV=0.64')),
'3_0.025_0.04'=
expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=2.5%,',italic('B'[u]),'=4% CV=0.63')),
'3_0.2_0.04'=
expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=20%,',italic('B'[u]),'=4% CV=0.62')),
'3_0.15_0.04'=
expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=15%,',italic('B'[u]),'=4% CV=0.61'))))+
guides(fill=guide_legend(ncol=2))+
geom_line(data=Yacorns, aes(y=(weight-mean(weight, na.rm = T))/sd(weight, na.rm = T)),
group=1)+
geom_point(data=Yacorns, aes(y=(weight-mean(weight, na.rm = T))/sd(weight, na.rm = T)))+
coord_cartesian(ylim=(c(-2.7, 5)))
## Warning: Removed 6 rows containing missing values (geom_point).
#RB+Pollen
Polim_kappa <-
Polim_f %>% filter(Year>1983) %>% select(parameter,Rc, Bu, Li, Year, run, masting) %>%
left_join(Yacorns[,c(1,8,9)], by='Year', suffix=c('_sim','_obs')) %>%
select(parameter, Rc, Bu, Li,run, Year, masting_sim, mastingW) %>%
group_by(parameter,Rc, Bu, Li, run) %>% nest() %>% #by run and parameter
mutate(data2=map(data, drop_y)) %>% #dropping year from nested frames
mutate(tables=map(data2,table)) %>% #getting into good table format
mutate(confKap2=map(tables, fun_Kap2), estimate=map(confKap2, est_Kap),
p=map(confKap2, p_Kap)) %>% convert(num(estimate,p)) %>%
group_by(parameter,Rc, Bu, Li) %>% summarise(est=mean(estimate)) %>%
mutate(Li=Li*100, Bu=Bu*100) %>%
convert(fct(Bu))
Polim_kappa %>% filter(Bu %in% c('3','4','5')) %>%
mutate(Bu=fct_relevel(Bu, "3",'4','5')) %>%
ggplot(aes(x=jitter(Li, factor = 0.4), y=est, color=Bu)) +
geom_point()+
geom_smooth(se=F, size=0.5)+
labs(title='(c)',x=expression(paste(italic('L'[T]),' (%)')),y='Kappa')+
scale_color_manual(values=c('black','blue','red'))+
theme_bw()+
theme(legend.position = 'none')+
coord_cartesian(ylim=c(-0.3,0.35))
#RB+Pollen
Polim_group <- Polim_f %>% filter(Year>2004 &Year <2018 &Year!=2010) %>%
group_by(parameter, Rc, Bu, Li, Year) %>% summarise(AcornW=mean(AcornW)) %>% #mean from 5 runs
left_join(Yseed, by=c('Year'='year')) %>%
mutate(error=(AcornW*1000*2.19)-(weight*0.8)) %>%
group_by(parameter, Rc, Bu, Li) %>%
summarise(RMSE=sqrt(mean(error^2, na.rm = T)),
MAE=mean(abs(error), na.rm=T)) %>% ungroup() %>%
mutate(Li=Li*100, Bu=Bu*100) %>%
select(Rc:RMSE) %>% convert(fct(Bu))
Polim_group %>% filter(Bu %in% c('3','4','5')) %>%
ggplot(aes(x=jitter(Li, factor = 0.4), y=RMSE, color=Bu)) +
geom_point(size=1)+
geom_smooth(se=F, size=0.5)+
labs(title='(c)',x=expression(paste(italic('L'[T]),' (%)')), y='RMSE (kg)')+
scale_color_manual(values=c('black','blue','red'))+
theme_bw()+
theme(legend.position = 'none')+
coord_cartesian(ylim=c(40,140))
Appendix S4
types <- c("Number" = 'dashed', "Mass" = 'solid')
Yacorns %>% mutate(CV=AcornSD/AcornMean, CV_W=AcornSDW/AcornMeanW) %>%
tail(16) %>%
ggplot(aes(x=Year))+
geom_line(aes(y=CV, linetype ='Number'))+
geom_line(aes(y=CV_W, linetype='Mass'))+
theme_classic()+
labs(x='End of 20-year window')+
theme(panel.grid.major.y = element_line(),
legend.position = 'none',
axis.text = element_text(size=13),
axis.title = element_text(size=15))+
scale_linetype_manual(values = types)+
scale_x_continuous(breaks=c(2003,2006, 2009, 2012, 2015, 2018))
Appendix S5
## Variability
#RB+Pollen
Polim_f %>% filter(Year>1983) %>%
group_by(Rc, Li, Bu, run) %>% summarise(CV=sd(AcornW)/mean(AcornW)) %>%
group_by(Rc, Li, Bu) %>% summarise(meanCV=mean(CV), seCV=sd(CV)/sqrt(5)) %>%
arrange(desc(meanCV)) %>% head(5)
## `summarise()` has grouped output by 'Rc', 'Li', 'Bu'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Rc', 'Li'. You can override using the `.groups` argument.
## Simultaneity
#RB+Pollen
Polim_f %>% filter(Year>1983) %>% select(parameter, Year, run, masting) %>%
left_join(Yacorns[,c(1,8,9)], by='Year', suffix=c('_sim','_obs')) %>%
select(parameter, run, Year, masting_sim, mastingW) %>%
group_by(parameter, run) %>% nest() %>% #by run and parameter
mutate(data2=map(data, drop_y)) %>% #dropping year from nested frames
mutate(tables=map(data2,table)) %>% #getting into good table format
mutate(confKap2=map(tables, fun_Kap2), estimate=map(confKap2, est_Kap),
p=map(confKap2, p_Kap)) %>% convert(num(estimate,p)) %>%
group_by(parameter) %>% summarise(est=mean(estimate)) %>%
arrange(desc(est)) %>% head(5)
## Effect size
Yacorn_freq <- Yacorns %>% filter(!is.na(freq)) %>%
select(Year, weight, freqW)
#RB+Pollen
Polim_f %>% filter(Year>2002) %>% #end of 20 year running window
left_join(Yacorn_freq) %>%
group_by(parameter, run) %>% nest() %>%
mutate(wilTest=map(data, wil_test_eff)) %>%
unnest(wilTest) %>% group_by(parameter) %>% summarise(eff=mean(effsize)) %>%
arrange(eff) %>% head(5)
## Joining, by = "Year"
## RMSE
#RB+Pollen
Polim_f %>% filter(Year>2004 &Year <2018 &Year!=2010) %>%
group_by(parameter, Year) %>% summarise(AcornW=mean(AcornW)) %>% #mean from 5 runs
left_join(Yseed, by=c('Year'='year')) %>%
mutate(error=(AcornW*1000*2.19)-(weight*0.8)) %>%
group_by(parameter) %>%
summarise(RMSE=sqrt(mean(error^2, na.rm = T))) %>%
arrange(RMSE) %>% head(5)
## `summarise()` has grouped output by 'parameter'. You can override using the `.groups` argument.